OHIBC data exploration: Mariculture

knitr::opts_chunk$set(fig.width = 6, fig.height = 4, fig.path = 'Figs/',
                      echo = TRUE, message = FALSE, warning = FALSE)

library(raster)
library(rgdal)
library(sf)
library(tmap)


source('~/github/ohibc/src/R/common.R')  ### an OHIBC specific version of common.R

scenario <- 'v2017'
goal     <- 'mar'
dir_git  <- path.expand('~/github/ohibc')
dir_goal <- file.path(dir_git, 'prep', goal, scenario)
dir_spatial <- file.path(dir_git, 'prep/_spatial')

dir_data_bc  <- file.path(dir_M, 'git-annex/bcprep', '_raw_data')

# library(provRmd); prov_setup()

### set up proj4string options: BC Albers and WGS84
p4s_bcalb <- c('bcalb' = '+init=epsg:3005')
p4s_wgs84 <- c('wgs84' = '+init=epsg:4326')

1 Summary: OHIBC Mariculture

This script prepares layers (employment rates and median income) for Livelihoods goal in British Columbia’s coastal regions.

From Halpern et al. (2014) (OHI California Current):

Livelihood sub-goal: As was done in the global analysis, coastal livelihoods is measured by two equally weighted sub-components, the number of jobs (j), which is a proxy for livelihood quantity, and the median annual household wages (g), which is a proxy for job quality. For jobs and wages we used a no-net loss reference point.

For British Columbia, we do not currently have sector-specific unemployment and wage information. As such we will analyze Livelihoods according to the model:

\(x_{LIV} = \frac{j' + g'}{2}\)

\(j' = \frac{j_c / j_{ref}}{M_c / M_{ref}}\)

where M is each region’s employment rate (1 - unemployment) as a percent at current (c) and reference (ref) time periods, and:

\(g' = \frac{g_c / g_{ref}}{W_c / W_{ref}}\)

where W is each State’s average annual per capita wage at current (c) and reference (ref) time periods.

---

2 Data sources

• Aquaculture potential data: Gentry/Froehlich 2017
• MaPP Aquaculture SMZs
• DFO Aquaculture tenures

---

3 Methods

3.1 Production potential

3.1.1 Convert Phi-prime raster to growth time

• Resample the \(\Phi'\) rasters from Gentry/Froehlich 2017 to OHIBC extents, 1 km resolution, and BC Albers projection.
• Convert \(\Phi'\) values to T_b and T_f values.
  • Fish: \(log(T_F) = 7.68 - 5.82 log(\Phi')\)
  • Bivalves: \(log(T_B) = 2.99 - 1.66 \Phi'\)

phi_rasts <- c('tf_rast' = file.path(dir_goal, 'data_explore/t_f_raster_1000m.tif'),
               'tb_rast' = file.path(dir_goal, 'data_explore/t_b_raster_1000m.tif'),
               'tb_minus_sd_rast' = file.path(dir_goal, 'data_explore/t_b_minus_sd_raster_1000m.tif'))
dir_aq <- file.path(dir_data_bc, 'aquaculture')

reload <- FALSE

if(any(!file.exists(phi_rasts)) | reload) {
  
  ohibc_rast <- raster(file.path(dir_spatial, 'raster/ohibc_rgn_raster_1000m.tif'))
  
  phif_rast <- raster(file.path(dir_aq, 'FishPhiALLConstraints95LT2.tif')) %>%
    projectRaster(ohibc_rast)
  
  tf_rast <- exp(7.86 - log(phif_rast) * 5.82)
  
  phib_rast <- raster(file.path(dir_aq, 'BivalvePhiALLConstraints95LT1.tif')) %>%
    projectRaster(ohibc_rast)
  
  phib_sd_rast <- raster(file.path(dir_aq, 'Bivalve_spp_Phi_sd.tif')) %>%
    projectRaster(ohibc_rast)
  
  # x <- phib_rast - phib_sd_rast
  
  tb_rast <- exp(2.99 - phib_rast * 1.66)
  
  tb_minus_sd_rast <- exp(2.99 - (phib_rast - phib_sd_rast) * 1.66)
  
  writeRaster(tf_rast, phi_rasts['tf_rast'], overwrite = TRUE)
  writeRaster(tb_rast, phi_rasts['tb_rast'], overwrite = TRUE)
  writeRaster(tb_minus_sd_rast, phi_rasts['tb_minus_sd_rast'], overwrite = TRUE)
}

phi_rasts_uncl <- c('tf_rast' = file.path(dir_goal, 'data_explore/t_f_unclipped_1000m.tif'),
               'tb_rast' = file.path(dir_goal, 'data_explore/t_b_unclipped_1000m.tif'),
               'tb_minus_sd_rast' = file.path(dir_goal, 'data_explore/t_b_minus_sd_unclipped_1000m.tif'))
dir_aq <- file.path(dir_data_bc, 'aquaculture')

reload <- FALSE

if(any(!file.exists(phi_rasts_uncl)) | reload) {
  ohibc_rast <- raster(file.path(dir_spatial, 'raster/ohibc_rgn_raster_1000m.tif'))
  
  phif_rast_uncl <- raster(file.path(dir_aq, 'spp_Phi_mean.tif')) %>%
    projectRaster(ohibc_rast)
  
  values(phif_rast_uncl)[values(phif_rast_uncl) < 2.0] <- NA
  ### clip out unreasonably low Phi values
  
  tf_rast_uncl <- exp(7.86 - log(phif_rast_uncl) * 5.82)
  
  phib_rast_uncl <- raster(file.path(dir_aq, 'Bivalve_spp_Phi_mean.tif')) %>%
    projectRaster(ohibc_rast)
  
  phib_sd_rast <- raster(file.path(dir_aq, 'Bivalve_spp_Phi_sd.tif')) %>%
    projectRaster(ohibc_rast)
  
  tb_rast_uncl <- exp(2.99 - phib_rast_uncl * 1.66)
  tb_minus_sd_rast_uncl <- exp(2.99 - (phib_rast_uncl - phib_sd_rast) * 1.66)
  
  writeRaster(tf_rast_uncl, phi_rasts_uncl['tf_rast'], overwrite = TRUE)
  writeRaster(tb_rast_uncl, phi_rasts_uncl['tb_rast'], overwrite = TRUE)
  writeRaster(tb_minus_sd_rast_uncl, phi_rasts_uncl['tb_minus_sd_rast'], overwrite = TRUE)
}

tf_rast <- raster(file.path(dir_goal, 'data_explore/t_f_raster_1000m.tif'))
tb_rast <- raster(file.path(dir_goal, 'data_explore/t_b_raster_1000m.tif'))
tf_rast_uncl <- raster(file.path(dir_goal, 'data_explore/t_f_unclipped_1000m.tif'))
tb_rast_uncl <- raster(file.path(dir_goal, 'data_explore/t_b_unclipped_1000m.tif'))

ohibc_sf <- sf::read_sf(dir_spatial, 'ohibc_rgn')
bc_land_sf <- sf::read_sf(dir_spatial, 'ohibc_land')

tf_map <- tm_shape(bc_land_sf) +
    tm_fill(col = 'grey40', alpha = 1) +
  tm_shape(ohibc_sf) +
    tm_fill(col = 'grey96', alpha = 1) +
  tm_shape(tf_rast_uncl) +
    tm_raster(title = 'Growth time\nyrs (35 cm fish)\n(all data)', 
              palette = 'Greens') +
  tm_shape(tf_rast) +
    tm_raster(title = 'Growth time\nyears (35 cm fish)', 
              palette = 'Reds') +
  tm_shape(ohibc_sf) +
    tm_borders(col = 'grey40', lwd = .25) +
  tm_legend(bg.alpha = .7, bg.color = 'white', position = c('right', 'top'))

print(tf_map)

tb_map <- tm_shape(bc_land_sf) +
    tm_fill(col = 'grey40', alpha = 1) +
  tm_shape(ohibc_sf) +
    tm_polygons(col = 'grey96', alpha = 1) +
  tm_shape(tb_rast_uncl) +
    tm_raster(title = 'Growth time\nyrs (4 cm bivalve)\n(all data)', 
              palette = 'Greens') +
  tm_shape(tb_rast) +
    tm_raster(title = 'Growth time\nyears (4 cm bivalve)',
              palette = 'Blues') +
  tm_shape(ohibc_sf) +
    tm_borders(col = 'grey40', lwd = .25) +
  tm_legend(bg.alpha = .7, bg.color = 'white', position = c('right', 'top'))

print(tb_map)

3.1.2 Calculate production biomass

From the growth time, we calculate biomass production based on the following assumptions and calculations:

• For fish, we assume (from Gentry/Froehlich) a cage stocking density (at harvest) of 11 kg/m^3, 9000 m^3 per cage, and 24 cages per km^2.
  • This results in a harvest per area of 2376 tonnes per km^2.
  • The rate of harvest is 1 harvest every \(1/T_F\) years.
  • Harvest intensity = \(2376 / T_F\) tonnes per year.
• For bivalves, we assume (from Gentry/Froehlich) 100 long lines per km^2, each of which contains 13,000 feet of fuzzy rope, seeded with 100 bivalves per foot.
  • This results in a harvest per area of 130e6 bivalves (4 cm) per km^2.
  • The rate of harvest is 1 harvest every \(1/T_B\) years.
  • Harvest intensity = \(130 x 10^6 / T_B\) bivalves per year.

tf_rast <- raster(file.path(dir_goal, 'data_explore/t_f_raster_1000m.tif'))
harvest_f_rast <- 2376/tf_rast

tb_rast <- raster(file.path(dir_goal, 'data_explore/t_b_raster_1000m.tif'))
harvest_b_rast <- 130e6/tb_rast

tb_sd_rast <- raster(file.path(dir_goal, 'data_explore/t_b_minus_sd_raster_1000m.tif'))
harvest_b_sd_rast <- 130e6/tb_sd_rast

writeRaster(harvest_f_rast, file.path(dir_goal, 'data_explore/harvest_f_raster_1000m.tif'), overwrite = TRUE)
writeRaster(harvest_b_rast, file.path(dir_goal, 'data_explore/harvest_b_raster_1000m.tif'), overwrite = TRUE)
writeRaster(harvest_b_sd_rast, file.path(dir_goal, 'data_explore/harvest_b_sd_raster_1000m.tif'), overwrite = TRUE)

tf_rast_uncl <- raster(file.path(dir_goal, 'data_explore/t_f_unclipped_1000m.tif'))
harvest_f_rast_uncl <- 2376/tf_rast_uncl

tb_rast_uncl <- raster(file.path(dir_goal, 'data_explore/t_b_unclipped_1000m.tif'))
harvest_b_rast_uncl <- 130e6/tb_rast_uncl

tb_sd_rast_uncl <- raster(file.path(dir_goal, 'data_explore/t_b_minus_sd_unclipped_1000m.tif'))
harvest_b_sd_rast_uncl <- 130e6/tb_sd_rast_uncl

writeRaster(harvest_f_rast_uncl, file.path(dir_goal, 'data_explore/harvest_f_unclipped_1000m.tif'), overwrite = TRUE)
writeRaster(harvest_b_rast_uncl, file.path(dir_goal, 'data_explore/harvest_b_unclipped_1000m.tif'), overwrite = TRUE)
writeRaster(harvest_b_sd_rast_uncl, file.path(dir_goal, 'data_explore/harvest_b_sd_unclipped_1000m.tif'), overwrite = TRUE)

harvest_f_rast <- raster(file.path(dir_goal, 'data_explore/harvest_f_raster_1000m.tif'))
harvest_b_rast <- raster(file.path(dir_goal, 'data_explore/harvest_b_raster_1000m.tif'))
harvest_b_sd_rast <- raster(file.path(dir_goal, 'data_explore/harvest_b_sd_raster_1000m.tif'))
harvest_f_rast_uncl <- raster(file.path(dir_goal, 'data_explore/harvest_f_unclipped_1000m.tif'), overwrite = TRUE)
harvest_b_rast_uncl <- raster(file.path(dir_goal, 'data_explore/harvest_b_unclipped_1000m.tif'), overwrite = TRUE)
harvest_b_sd_rast_uncl <- raster(file.path(dir_goal, 'data_explore/harvest_b_sd_unclipped_1000m.tif'), overwrite = TRUE)

ohibc_sf <- sf::read_sf(dir_spatial, 'ohibc_rgn')
bc_land_sf <- sf::read_sf(dir_spatial, 'ohibc_land')

harvest_f_map <- tm_shape(bc_land_sf) +
    tm_fill(col = 'grey40', alpha = 1) +
  tm_shape(ohibc_sf) +
    tm_fill(col = 'grey80', alpha = .4) +
  tm_shape(harvest_f_rast_uncl) +
    tm_raster(title = 'Fish harvest\ntonnes/yr (all pts)', 
              palette = 'Greens') +
  tm_shape(harvest_f_rast) +
    tm_raster(title = 'Fish harvest\ntonnes/year') +
  tm_shape(ohibc_sf) +
    tm_borders(col = 'grey40', lwd = .25) +
  tm_legend(bg.alpha = .7, bg.color = 'white', position = c('right', 'top'))

print(harvest_f_map)

harvest_b_map <- tm_shape(bc_land_sf) +
    tm_fill(col = 'grey40', alpha = 1) +
  tm_shape(ohibc_sf) +
    tm_polygons(col = 'grey80', alpha = .4) +
  tm_shape(harvest_b_rast_uncl) +
    tm_raster(title = 'Bivalve harvest\nunits/yr (all pts)', 
              palette = 'Greens') +
  tm_shape(harvest_b_rast) +
    tm_raster(title = 'Bivalve harvest\nunits/yr',
              palette = 'Blues') +
  tm_shape(ohibc_sf) +
    tm_borders(col = 'grey40', lwd = .25) +
  tm_legend(bg.alpha = .7, bg.color = 'white', position = c('right', 'top'))

print(harvest_b_map)

harvest_b_sd_map <- tm_shape(bc_land_sf) +
    tm_fill(col = 'grey40', alpha = 1) +
  tm_shape(ohibc_sf) +
    tm_polygons(col = 'grey80', alpha = .4) +
  tm_shape(harvest_b_sd_rast_uncl) +
    tm_raster(title = 'Bivalve harvest\nunits/yr (-1 sd)', 
              palette = 'Greens') +
  tm_shape(harvest_b_sd_rast) +
    tm_raster(title = 'Bivalve harvest\nunits/yr (-1 sd)',
              palette = 'Blues') +
  tm_shape(ohibc_sf) +
    tm_borders(col = 'grey40', lwd = .25) +
  tm_legend(bg.alpha = .7, bg.color = 'white', position = c('right', 'top'))

print(harvest_b_sd_map)

3.1.3 Explore data by region

For each OHIBC region, examine the distribution of production potential for both fish and bivalves. Some possible targets could be based on median production for each region, other quantiles. Due to coarseness of production raster, production hot spots do not line up well with tenure locations, so spatially estimating production potential based on tenures is not likely to be a good method.

• For finfish aquaculture, which occurs primarily in the southern regions, sufficient data points exist in the regions to determine a reasonable cross-section of representative values.
  
• For shellfish aquaculture, the clipped \(Phi'\) dataset cuts all cells, probably due to insufficient chlorophyll to supply the bivalves according to the chlorophyll data. However, we will use the unclipped \(Phi'\) data to identify a range of reasonable production values, assuming that producers figure out how to identify locations with sufficient chlorophyll.
• Since shellfish aquaculture potential seems far beyond actual production, we will also explore a lower bound to account for uncertainty in shellfish growth calculations. For this, we will take \(\Phi' - sd(\Phi')\) for each cell.

harvest_f_rast <- raster(file.path(dir_goal, 'data_explore/harvest_f_raster_1000m.tif'))

### use extract() to get values per region
ohibc_poly <- readOGR(dir_spatial, 'ohibc_rgn')

## OGR data source with driver: ESRI Shapefile 
## Source: "/home/ohara/github/ohibc/prep/_spatial", layer: "ohibc_rgn"
## with 8 features
## It has 5 fields
## Integer64 fields read as strings:  rgn_id

f_values <- raster::extract(harvest_f_rast, ohibc_poly) %>%
  lapply(FUN = function(x) data.frame(prod_potential = x)) %>%
  setNames(ohibc_poly$rgn_id) %>%
  bind_rows(.id = 'rgn_id') %>%
  mutate(rgn_id = as.integer(rgn_id),
         prod_potential = round(prod_potential, 3))

prob_levels <- c(.01, .02, .05, .10, .25, .50, .75, .90, .95, .98, .99)
f_range <- f_values %>%
  group_by(rgn_id) %>%
  summarize(n_tot = n(), n_phi = sum(!is.na(prod_potential)), 
            min_prod = min(prod_potential, na.rm = TRUE), max_prod = max(prod_potential, na.rm = TRUE),
            quants = list(stats::quantile(prod_potential, probs = prob_levels, na.rm = TRUE))) %>%
  unnest() %>%
  mutate(prob_lvl = rep(prob_levels, times = 8)) %>%
  ungroup()

f_plot <- ggplot(f_values %>%
                   left_join(get_rgn_names(), by = 'rgn_id')) + 
  ggtheme_plot() +
  geom_histogram(aes(x = prod_potential)) + 
  facet_wrap(~ rgn_name, scales = 'free_y') +
  labs(x = 'Production potential (Finfish tonnes/km^2/year)',
       y = 'Number of 1 km^2 cells',
       title = 'Finfish')
print(f_plot)

write_csv(f_range, file.path(dir_goal, 'data_explore/prod_pot_f_range.csv'))

harvest_b_rast <- raster(file.path(dir_goal, 'data_explore/harvest_b_raster_1000m.tif'))

b_values <- raster::extract(harvest_b_rast, ohibc_poly) %>%
  lapply(FUN = function(x) data.frame(prod_potential = x)) %>%
  setNames(ohibc_poly$rgn_id) %>%
  bind_rows(.id = 'rgn_id') %>%
  mutate(rgn_id = as.integer(rgn_id),
         prod_potential = round(prod_potential))

prob_levels <- c(.01, .02, .05, .10, .25, .50, .75, .90, .95, .98, .99)
b_range <- b_values %>%
  group_by(rgn_id) %>%
  summarize(n_tot = n(), n_phi = sum(!is.na(prod_potential)), 
            min_prod = min(prod_potential, na.rm = TRUE), max_prod = max(prod_potential, na.rm = TRUE),
            quants = list(stats::quantile(prod_potential, probs = prob_levels, na.rm = TRUE))) %>%
  unnest() %>%
  mutate(prob_lvl = rep(prob_levels, times = 8)) %>%
  ungroup()

b_plot <- ggplot(b_values %>%
                   left_join(get_rgn_names(), by = 'rgn_id')) + 
  ggtheme_plot() +
  geom_histogram(aes(x = prod_potential)) + 
  facet_wrap(~ rgn_name, scales = 'free_y') +
  labs(x = 'Production potential (4 cm bivalves/km^2/year)',
       y = 'Number of 1 km^2 cells',
       title = 'Bivalves')
print(b_plot)

write_csv(b_range, file.path(dir_goal, 'data_explore/prod_pot_b_range.csv'))

In generating median harvest values for bivalves from the unclipped dataset (to fill the NAs in North Coast, and flesh out a few more datapoints for Strait of Georgia and Haida Gwaii), we want to ensure that the numbers aren’t outrageous. For the regions with data in the clipped dataset, we compare to the unclipped and see that the medians still fall fairly close.

harvest_b_uncl_rast <- raster(file.path(dir_goal, 'data_explore/harvest_b_unclipped_1000m.tif'))

b_uncl_values <- raster::extract(harvest_b_uncl_rast, ohibc_poly) %>%
  lapply(FUN = function(x) data.frame(prod_potential = x)) %>%
  setNames(ohibc_poly$rgn_id) %>%
  bind_rows(.id = 'rgn_id') %>%
  mutate(rgn_id = as.integer(rgn_id),
         prod_potential = round(prod_potential))

prob_levels <- c(.01, .02, .05, .10, .25, .50, .75, .90, .95, .98, .99)
b_uncl_range <- b_uncl_values %>%
  group_by(rgn_id) %>%
  summarize(n_tot = n(), n_phi = sum(!is.na(prod_potential)), 
            min_prod = min(prod_potential, na.rm = TRUE), max_prod = max(prod_potential, na.rm = TRUE),
            quants = list(stats::quantile(prod_potential, probs = prob_levels, na.rm = TRUE))) %>%
  unnest() %>%
  mutate(prob_lvl = rep(prob_levels, times = 8)) %>%
  ungroup()

b_uncl_plot <- ggplot(b_uncl_values %>%
                   left_join(get_rgn_names(), by = 'rgn_id')) + 
  ggtheme_plot() +
  geom_histogram(aes(x = prod_potential)) + 
  facet_wrap(~ rgn_name, scales = 'free_y') +
  labs(x = 'Production potential (4 cm bivalves/km^2/year) (all data)',
       y = 'Number of 1 km^2 cells',
       title = 'Bivalves (all data)')
print(b_uncl_plot)

write_csv(b_uncl_range, file.path(dir_goal, 'data_explore/prod_pot_b_range_unclipped.csv'))

x <- b_range %>%
  select(rgn_id, quants, prob_lvl) %>%
  left_join(b_uncl_range %>%
              select(rgn_id, quants_unclipped = quants, prob_lvl),
            by = c('rgn_id', 'prob_lvl')) %>%
  # filter(prob_lvl == .50) %>%
  left_join(get_rgn_names())

test_plot <- ggplot(x, aes(x = quants, y = quants_unclipped, label = prob_lvl)) +
  geom_abline(intercept = 0, slope = 1, color = 'red') +
  geom_hline(data = x %>% filter(is.na(quants) & prob_lvl == 0.5), 
             aes(yintercept = quants_unclipped, label = rgn_name), alpha = .5) +
  geom_point(alpha = .3, show.legend = FALSE) +
  geom_point(data = x %>% filter(prob_lvl == 0.5), aes(color = rgn_name)) +
  xlim(0, NA) + ylim(0, NA) +
  labs(title = 'Median bivalve production, \nall data vs clipped',
       x = 'prod. potential clipped',
       y = 'prod. potential, all cells')

plotly::ggplotly(test_plot)

Since there is such uncertainty around bivalve growth rates let’s examine the impact of subtracting the standard deviation, as a possible lower target. This plot considers unclipped values for all regions.

harvest_b_sd_uncl_rast <- raster(file.path(dir_goal, 'data_explore/harvest_b_sd_unclipped_1000m.tif'))

b_sd_uncl_values <- raster::extract(harvest_b_sd_uncl_rast, ohibc_poly) %>%
  lapply(FUN = function(x) data.frame(prod_potential = x)) %>%
  setNames(ohibc_poly$rgn_id) %>%
  bind_rows(.id = 'rgn_id') %>%
  mutate(rgn_id = as.integer(rgn_id),
         prod_potential = round(prod_potential))

prob_levels <- c(.01, .02, .05, .10, .25, .50, .75, .90, .95, .98, .99)
b_sd_uncl_range <- b_sd_uncl_values %>%
  group_by(rgn_id) %>%
  summarize(n_tot = n(), n_phi = sum(!is.na(prod_potential)), 
            min_prod = min(prod_potential, na.rm = TRUE), max_prod = max(prod_potential, na.rm = TRUE),
            quants = list(stats::quantile(prod_potential, probs = prob_levels, na.rm = TRUE))) %>%
  unnest() %>%
  mutate(prob_lvl = rep(prob_levels, times = 8)) %>%
  ungroup()

b_sd_uncl_plot <- ggplot(b_sd_uncl_values %>%
                   left_join(get_rgn_names(), by = 'rgn_id')) + 
  ggtheme_plot() +
  geom_histogram(aes(x = prod_potential)) + 
  facet_wrap(~ rgn_name, scales = 'free_y') +
  labs(x = 'Production potential (4 cm bivalves/km^2/year) (-1 sd)',
       y = 'Number of 1 km^2 cells',
       title = "Bivalves (median based on Phi' - 1sd)")
print(b_sd_uncl_plot)

write_csv(b_sd_uncl_range, file.path(dir_goal, 'data_explore/prod_pot_b_sd_range_unclipped.csv'))

x <- b_uncl_range %>%
  select(rgn_id, quants, prob_lvl) %>%
  left_join(b_sd_uncl_range %>%
              select(rgn_id, quants_minus_sd = quants, prob_lvl),
            by = c('rgn_id', 'prob_lvl')) %>%
  # filter(prob_lvl == .50) %>%
  left_join(get_rgn_names())

test_plot <- ggplot(x, aes(x = quants, y = quants_minus_sd, label = prob_lvl)) +
  geom_abline(intercept = 0, slope = 1, color = 'red') +
  geom_point(alpha = .3, show.legend = FALSE) +
  geom_point(data = x %>% filter(prob_lvl == 0.5), aes(color = rgn_name)) +
  xlim(0, NA) + ylim(0, NA) +
  labs(title = 'Bivalve production, -1 sd',
       x = 'prod. potential, all cells',
       y = "prod. potent, phi' - 1sd")

plotly::ggplotly(test_plot)

3.2 Aquaculture development area

As aquaculture development targets, we will use MaPP Special Management Zones assigned to aquaculture to represent the goals of the MaPP regions. For non-MaPP regions, we will use currently allocated aquaculture tenures (active/not, developed/not) to indicate development goals.

3.2.1 MaPP Special Management Zones for aquaculture

To determine the area allocated to aquaculture/mariculture for the MaPP regions, we will use SMZ proportions of total region ocean area, taken from MaPP region shapefiles. In some cases (e.g. Haida Gwaii west coast north, table 8.21) it appears some non-SMZ areas are noted “acceptable” for aquaculture. As a conservative estimate, only SMZ areas will be counted.

Note that for North Coast Vancouver Island, the OHIBC region wraps around to the west in addition to the overall MaPP region. In this case, perhaps we want to include existing aquaculture tenure areas? For now, these are excluded.

Do any of these regions allow finfish aquaculture? If so how to allocate the area appropriately…

• Haida Gwaii and NC MaPP: shellfish only (no production in NC 2011-2015); CC not specified but currently (2011-2015) only produces finfish; and NCVI states (and produces) both, but MaPP regions are explicitly for shellfish.

dir_mapp <- file.path(dir_data_bc, 'mapp/MarineSpatialPlanZones')

mapp_hg   <- read_sf(dir_mapp, 'HaidaGwaii_Subregion_Oct_2014')
    # "SubRegion"  "Type"       "Management" "Area_Ha"    "Area_Km"    "Objective"  "Name"       
    # "Haida_Name" "Shape_Leng" "Shape_Area"
    # mapp_hg$Management %>% unique()
    # [1] "SMZ - Shellfish Aq" "IUCN - Type Ib"     "IUCN - Type IV"     "IUCN - Type II"    
    # [5] "IUCN - Type III"    "IUCN - Type V"      "IUCN - Type VI"     "SMZ - Alt. Energy" 
    mapp_hg_aq <- mapp_hg %>%
      filter(Management == 'SMZ - Shellfish Aq') %>%
      mutate(area = st_area(geometry)) %>%
      as.data.frame() %>%
      dplyr::select(name = Name, area, desc = Management) %>%
      mutate(rgn_id = 2)
mapp_nc   <- read_sf(dir_mapp, 'NorthCoast_Subregion_June2015')
    # "SubRegion"  "Type"       "Management" "Area_Ha"    "Area_Km"    "Objective"  "AreaDescri" 
    # "Name"       "Unit_No"    "Grouping"   "ZoneType"   "Shape_Leng" "Shape_Area"
    # mapp_nc$Management %>% unique()
    # [1] "II"                                      "IV"                                     
    # [3] "Renewable Energy"                        "Aquaculture"                            
    # [5] "Tourism and Recreation"                  "Areas for future planning consideration"
    # [7] "Ib"                                      "Cultural" 
    mapp_nc_aq <- mapp_nc %>%
      filter(Management == 'Aquaculture') %>%
      mutate(area = st_area(geometry)) %>%
      as.data.frame() %>%
      dplyr::select(name = Name, area, desc = Management) %>%
      mutate(rgn_id = 1)
mapp_cc   <- read_sf(dir_mapp, 'CentralCoast_Subregion_Ver9')
    # "SubRegion"  "Type"       "Area_Ha"    "Area_Km"    "Unit_New"   "Unit"       "Group_"    
    # "IUCN_Categ" "ZoneType"   "Edited"     "Unit_Old"   "AreaDecri"  "Shape_Leng" "Shape_Area"
    # mapp_cc$ZoneType %>% unique()
    # [1] NA                       "Renewable Energy"       "Aquaculture"           
    # [4] "Recreation and Tourism"
    mapp_cc_aq <- mapp_cc %>%
      filter(ZoneType == 'Aquaculture') %>%
      mutate(area = st_area(geometry)) %>%
      as.data.frame() %>%
      dplyr::select(name = AreaDecri, area, desc = ZoneType) %>%
      mutate(rgn_id = 3)
mapp_ncvi <- read_sf(dir_mapp, 'NorthVancouverIsland_Subregion_Oct3_2014')
    # "SubRegion"  "Type"       "Management" "Area_Ha"    "Area_Km"    "AreaDescri" "Name"      
    # "Unit_No"    "Shape_Leng" "Shape_Area"
    # mapp_ncvi$AreaDescri %>% unique()
    mapp_ncvi_aq <- mapp_ncvi %>%
      filter(str_detect(tolower(AreaDescri), 'aquaculture')) %>%
      mutate(area = st_area(geometry)) %>%
      as.data.frame() %>%
      dplyr::select(name = Name, area, desc = AreaDescri) %>%
      mutate(rgn_id = 4)
    # [1] "SMZ Cultural/Economic Emphasis Areas are intended to reinforce their high value to First Nations, on a seasonal and year-round basis, for cultural value protection, Aboriginal economic development opportunities, and food security. This value includes con"
    #  [2] "Significant ecological values due to major upwelling of nutrients creating a rich, diverse marine ecosystem.  There are key First Nation cultural/economic interests and local resident scenic values. Safeguarding the integrity of this interaction between" 
    #  [3] "SMZ Recreation/Tourism Emphasis Areas are intended to reinforce their high value to existing commercial tourism operations, particularly during the months of late May to early October. Other uses and activities in SMZ Recreation/Tourism Emphasis Areas sh"
    #  [4] "SMZ Community Emphasis Areas are intended to reinforce their value for seasonal and year-round uses and activities associated with, required by, and primarily dictated by, adjacent, or nearby communities. The uses and activities in SMZ Community Emphasis"
    #  [5] "SMZ Shellfish Aquaculture Emphasis Areas are intended to reinforce interest by First Nations in investigation and (if feasible) the development of bottom and off-bottom shellfish aquaculture operations. These areas may be associated with integrated multi"
    #  [6] "Important species and habitats, including those of cultural importance to First Nations. Significant for seasonal marine harvesting and ecotourism activities by First Nations. It is an important whale and wildlife viewing area. Includes important habitat"
    #  [7] "Important habitat and species, in particular a significant and unique glass sponge reef formation, which includes a complex ecosystem, enabling a species-rich marine environment that supports the local biodiversity of the area. Contains critical habitat" 
    #  [8] "The area is representative of shallow sill ecosystems of coral fans, sponges. Several deepwater and/or rare species including the gorgonian coral, the soft goblet sponge, the cloud sponge, the townsend eualid shrimp, and the bigmouth sculpin are found at"
    #  [9] "Considerable cultural modification by First Nations based on use of important marine species and habitat.  Ongoing practices and teachings, restoration of First Nations’ cultural resources, and their associated marine species and habitats, and for repa"
    # [10] "Important marine species, habitats, First Nations’ cultural resources such as cultural tourism, loxiwe, shell middens, and former seasonal village/resource processing site."
    # [11] "High marine recreational values, containing important marine species and habitats including important areas for herring and northern resident killer whales. Area includes First Nations’ cultural resources uses and activities such as cultural tourism, l"  
    # [12] "A diverse marine ecosystem, with important marine species and habitat. Important recreation and tourism area which includes several scuba diving sites. Includes important areas for humpback and and northern resident killer whales, herring and sea otters."
    # [13] "Marine species and habitats including those of cultural importance to First Nations. Connects existing conservation and protection areas and provides network/corridor between the Central Coast and NVI marine plans to assist in conservation and protection"
    # [14] "Important marine species and habitats including herring important areas. Protection of representative marine ecosystems at the confluence of three channels supporting rich intertidal species and habitats."

mapp_aq_area_df <- bind_rows(mapp_hg_aq, 
                             mapp_cc_aq, 
                             mapp_nc_aq, 
                             mapp_ncvi_aq) %>%
  mutate(area_km2 = round(area / 1e6, 3)) %>%
  select(-area)

write_csv(mapp_aq_area_df, file.path(dir_goal, 'data_explore/mapp_aq_smz_areas.csv'))

3.2.2 DFO Aquaculture tenures

For Strait of Georgia and West Coast Vancouver Island, we will use the total area of designated tenures according to DFO records, split into finfish and shellfish.

Note that for Aristazabal Island OHIBC region, no aquaculture areas (tenures or SMZs) have been identified.

tenure_shps <- list.files(file.path(dir_data_bc, 'dfo_khunter/aquaculture/d2016/TENURES_PNT'),
                          pattern = '.shp$', full.names = TRUE)
ohibc_rgn <- read_sf(dir_spatial, 'ohibc_rgns_unclipped')

tenures_ohibc <- lapply(tenure_shps, FUN = function(x) {
  ### x <- tenure_shps[1]
  read_sf(dsn = dirname(x), layer = str_replace(basename(x), '.shp$', '')) %>%
    mutate(area = st_area(geometry)) %>%
    st_intersection(ohibc_rgn) %>%
    as.data.frame() %>%
    dplyr::select(-geometry)
}) %>%
  setNames(str_replace(basename(tenure_shps), '.shp$', '')) %>%
  bind_rows(.id = 'filename')
  
tenure_areas <- tenures_ohibc %>%
  setNames(tolower(names(.))) %>% 
  mutate(spp = ifelse(is.na(sp_lic), sp_, sp_lic),
         area_km2 = round(area / 1e6, 3)) %>%
  dplyr::select(filename, rgn_id, location, spp, area_km2)

write_csv(tenure_areas, file.path(dir_goal, 'data_explore/dfo_aq_tenure_areas.csv'))

3.3 Get regional aquaculture harvests

DFO data for aquaculture harvests by PFMA for 2011-2015. Harvest values are in tonnes?

aq_dfo_file <- file.path(dir_data_bc, 'dfo_khunter/aquaculture/d2016',
                         'Aquaculture_production_PFMA 2011-15.xlsx')
data_aq_dfo <- readxl::read_excel(aq_dfo_file, skip = 4) %>%
  setNames(tolower(names(.)) %>%
             str_replace_all('[^a-z]+', '_'))

### This file has a number of rows for shellfish, then a number of rows for
### finfish.  The header row can be detected by looking for non-numeric
### text in the "year" column.
start_fish <- which(data_aq_dfo$year %>% str_detect('[a-zA-Z]'))

aq_shellfish <- data_aq_dfo[1:(min(start_fish) - 1), ] %>%
  gather(key = 'species', value = 'harvest', -year, -pfma) %>%
  mutate(aq_type = 'shellfish')

finfish_hdr  <- as.character(data_aq_dfo[max(start_fish), ]) %>%
  tolower() %>%
  str_replace_all('[^a-z]+', '_') %>%
  .[. != 'na']

aq_finfish   <- data_aq_dfo[(max(start_fish) + 1):nrow(data_aq_dfo), 1:length(finfish_hdr)] %>%
  setNames(finfish_hdr) %>%
  gather(key = 'species', value = 'harvest', -year, -pfma) %>%
  mutate(aq_type = 'finfish')

aq_harvest_df <- bind_rows(aq_shellfish, aq_finfish) %>%
  filter(!is.na(year))
  

write_csv(aq_harvest_df, file.path(dir_goal, 'data_explore/dfo_aq_harvest_pfma.csv'))

Since the DFO harvest is provided by PFMA, collapse PFMAs to OHIBC regions. We have already done this in the FIS goal, so we will use the same table to assign PFMAs here.

pfma_to_ohibc <- read_csv(file.path(dir_git, 'prep/fis/v2017/int/pfmsa_to_ohibc.csv')) %>%
  group_by(pfma_id, rgn_id) %>%
  summarize(area_km2 = sum(area_km2)) %>%
  group_by(pfma_id) %>%
  mutate(prop_area = area_km2 / sum(area_km2)) %>%
  select(pfma = pfma_id, rgn_id, prop_area)
  
aq_harvest_ohibc <- read_csv(file.path(dir_goal, 'data_explore/dfo_aq_harvest_pfma.csv')) %>%
  left_join(pfma_to_ohibc, by = 'pfma') %>%
  mutate(harvest_rgn = harvest * prop_area) %>%
  select(year, species, harvest_rgn, aq_type, rgn_id) %>%
  group_by(year, species, aq_type, rgn_id) %>%
  summarize(harvest = sum(harvest_rgn)) %>%
  ungroup()

write_csv(aq_harvest_ohibc, file.path(dir_goal, 'data_explore/dfo_aq_harvest_ohibc.csv'))

DT::datatable(aq_harvest_ohibc)

aq_plot_df <- aq_harvest_ohibc %>%
         left_join(get_rgn_names(), by = 'rgn_id') %>%
  mutate(species = str_replace(species, '_', ' '),
         species = str_replace(species, '_[a-z_]+', ''),
         species = ifelse(str_detect(species, 'sable'), 'sablefish', species))

aq_sum_df <- aq_plot_df %>%
  group_by(year, rgn_name, aq_type) %>%
  summarize(harvest = sum(harvest)) %>%
  ungroup()

ggplot(aq_plot_df,
       aes(x = year, y = harvest, color = species)) +
  ggtheme_plot() +
  geom_line(data = aq_sum_df %>% filter(aq_type == 'finfish'),
            aes(x = year, y = harvest), 
            size = 1.5, linetype = 'dashed', color = 'grey30') +
  geom_line(data = aq_sum_df %>% filter(aq_type == 'shellfish'),
            aes(x = year, y = harvest), 
            size = 1.5, linetype = 'dotted', color = 'grey30') +
  geom_point() +
  geom_line(aes(group = species)) +
  facet_wrap( ~ rgn_name, scale = 'free_y')

3.4 Create and write layers

Layers to ~/github/ohibc/prep/mar/v2017/output:

• Aquaculture development area by region (finfish and shellfish)
  • Note: CC MaPP area is allocated to finfish, as only finfish have current (2011-2015) production. All other MaPP areas are allocated to shellfish by explicit description.
• Aquaculture median productivity by region (finfish and shellfish)
• Aquaculture harvest by region (finfish and shellfish)

mapp_aq_types <- data.frame(rgn_id  = c(1:4),
                            aq_type = c('shellfish', 'shellfish', 'finfish', 'shellfish'))

mapp_aq_areas <- read_csv(file.path(dir_goal, 'data_explore/mapp_aq_smz_areas.csv')) %>%
  group_by(rgn_id, desc) %>%
  summarize(area_km2 = sum(area_km2)) %>%
  ungroup() %>%
  left_join(mapp_aq_types, by = 'rgn_id') %>%
  select(rgn_id, aq_type, area_km2) %>%
  mutate(source = 'mapp')
  
dfo_aq_areas <- read_csv(file.path(dir_goal, 'data_explore/dfo_aq_tenure_areas.csv')) %>%
  mutate(aq_type = ifelse(str_detect(filename, 'Finfish'), 'finfish', 'shellfish')) %>%
  group_by(rgn_id, aq_type) %>%
  summarize(area_km2 = sum(area_km2)) %>%
  ungroup() %>%
  distinct() %>%
  mutate(source = 'dfo')


aq_areas <- bind_rows(mapp_aq_areas, dfo_aq_areas) %>%
  group_by(rgn_id, aq_type) %>%
  mutate(a_tot_km2 = sum(area_km2)) %>%
  ungroup()
  
write_csv(aq_areas, file.path(dir_goal, 'output/aq_areas.csv'))

DT::datatable(aq_areas)

To compare potential vs harvest, we need to convert bivalve units to metric tonnes. Some figures:

• Scallops: 40-50 per kg (~22g ea)
• Mussels: 15-20 per 1 lb blue (~25g ea)
• Oysters: 12-15 per 1 lb (~33g ea) aquaculture halfshell
• Clams: 20-30 per 1 lb Manila (~18g ea), 8-14 per 1 lb littlenecks (~40g ea)

Averaging these gives about 27.5g per piece

aq_mass_per_pc <- 0.0275e-3 ### tonnes per piece

pot_b_unclipped <- read_csv(file.path(dir_goal, 'data_explore/prod_pot_b_range_unclipped.csv')) %>%
  filter(prob_lvl == .50) %>%
  select(rgn_id, median_prod_uncl = quants)
pot_b <- read_csv(file.path(dir_goal, 'data_explore/prod_pot_b_range.csv')) %>%
  filter(prob_lvl == .50) %>%
  select(rgn_id, median_prod = quants) %>%
  left_join(pot_b_unclipped, by = 'rgn_id') %>%
  mutate(median_prod_tonnes = ifelse(is.na(median_prod), median_prod_uncl, median_prod),
         median_prod_tonnes = median_prod_tonnes * aq_mass_per_pc,
         aq_type = 'shellfish',
         units   = 'tonnes') %>%
  select(rgn_id, median_prod = median_prod_tonnes, aq_type, units)

pot_f <- read_csv(file.path(dir_goal, 'data_explore/prod_pot_f_range.csv')) %>%
  filter(prob_lvl == 0.50) %>%
  select(rgn_id, median_prod = quants) %>%
  mutate(aq_type = 'finfish',
         units   = 'tonnes')

pot_aq <- bind_rows(pot_b, pot_f)

write_csv(pot_aq, file.path(dir_goal, 'output/aq_potential.csv'))

DT::datatable(pot_aq)

harvest_df <- read_csv(file.path(dir_goal, 'data_explore/dfo_aq_harvest_ohibc.csv')) %>%
  group_by(year, rgn_id, aq_type) %>%
  summarize(harvest_tonnes = sum(harvest)) %>%
  ungroup()

write_csv(harvest_df, file.path(dir_goal, 'output/aq_harvest.csv'))

DT::datatable(harvest_df)

3.5 Visualize

Reference points pictured are based only on DFO tenures; the MaPP reference areas are far larger and would create far more ambitious reference points for shellfish.

aq_area      <- read_csv(file.path(dir_goal, 'output', 'aq_areas.csv'))
aq_potential <- read_csv(file.path(dir_goal, 'output', 'aq_potential.csv'))
aq_harvest   <- read_csv(file.path(dir_goal, 'output', 'aq_harvest.csv'))


aq_all <- aq_area %>%
  left_join(aq_potential, by = c('rgn_id', 'aq_type')) %>%
  left_join(aq_harvest, by = c('rgn_id', 'aq_type')) %>%
  mutate(ref_pt = area_km2 * median_prod) %>%
  left_join(get_rgn_names(), by = 'rgn_id')

aq_f_df <- aq_all %>%
  filter(aq_type == 'finfish') %>%
  filter(source == 'dfo') %>%
  select(year, rgn_name, harvest_tonnes, ref_pt)

aq_f_plot <- ggplot(aq_f_df, aes(x = year, y = harvest_tonnes)) +
  ggtheme_plot() +
  geom_line(aes(group = rgn_name)) +
  geom_hline(aes(yintercept = ref_pt), color = 'red', alpha = .5) +
  labs(title = 'Finfish harvest',
       color = 'Ref pt source') +
  facet_wrap( ~ rgn_name, scales = 'free_y') +
  ylim(0, NA)

print(aq_f_plot)

aq_b_df <- aq_all %>%
  filter(aq_type == 'shellfish') %>%
  filter(source == 'dfo') %>%
  select(year, rgn_name, harvest_tonnes, ref_pt)

aq_b_plot <- ggplot(aq_b_df, aes(x = year, y = harvest_tonnes)) +
  ggtheme_plot() +
  geom_line(aes(group = rgn_name)) +
  geom_hline(aes(yintercept = ref_pt), color = 'red', alpha = .5) +
  labs(title = 'Shellfish harvest',
       color = 'Ref pt source') +
  facet_wrap( ~ rgn_name, scales = 'free_y') +
  ylim(0, NA)

print(aq_b_plot)

---

# prov_wrapup()